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ABSTRACT 

Gamma-ray burst (GRB) afterglows arc well described by synchrotron emission origi- 
nating from the interaction between a relativistic blast wave and the external medium 
surrounding the GRB progenitor. We introduce a code to reconstruct spectra and light 
curves from arbitrary fluid configurations, making it especially suited to study the ef- 
fects of fluid flows beyond those that can be described using analytical approximations. 
As a check and first application of our code we use it to fit the scaling coefficients of 
theoretical models of afterglow spectra. We extend earlier results of other authors 
to general circumburst density profiles. We rederive the physical parameters of GRB 
970508 and compare with other authors. 
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1 INTRODUCTION 

In the fireball model, Gamma-Ray Burst (GRB) afterglows 
are thought to be the result of synchrotron radiation gener- 
ated by electrons during the interaction of a strongly colli- 
mated relativistic jet from a comp act source w i th its en- 
viron ment (for recent reviews, see IPiranI 120051 ; IMeszarosI 
2006). Initially the resulting spectra and light curves have 
been modelled using only the shock front of a spherical 
explosion and a simple po wer law approximation for the 
synchrotron radiation (e.g. IWiiers. Rees fc M cszaros 1 19971 : 
iMeszaros fc Reeslll997l ; ISari et al.lll998l : iRhoadsll 19991 ). One 
or more spectral and temporal breaks were used to connect 
regimes with different power law slopes. For the dynamics 
the se lf similar approximation of a relativistic explosion was 
used l|Blandford fc McKeel fl976l ) . These models have been 
refined continuou sly. More details of the shock structure 
were included (e.g. lGranot et al. 1999; Gruzinov fc Waxmanl 
1999), more accur ate formulae for the sy nchrotron radia- 
tion were used (e.g. lWiiers fc GalamallT999l ) and efforts have 
been made to implement collimation using various analyti- 
cal approxima tions to the je t structure and lateral spreading 
behaviour fsee lGranotil2005l for an overview). On top of that, 
there have been stu dies focussing on arrival time effects (e.g. 
Huang et all 120071) and some numerical simulations (e.g . 
Salmonsonll2003l : iGranot et aill200ll : iNakar fc Granotll2007V) . 

The aim of this paper is twofold. The first aim is to in- 
troduce a new method to derive light curves and spectra by 
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post-processing relativistic hydrodynamic (RHD) jet simu- 
lations of arbitrary dimension, properly taking into account 
all beaming and arrival time effects, as well as the precise 
shape of the synchrotron spectrum and electron cooling (in 
this paper we will ignore self-absorption, although it can in 
principle be included in our method). This is done in sections 
II and O 

The second aim is to present a set of scaling coeffi- 
cients for the slow-cooling case for a density profile p — po ■ 
(R/Ro)~ k for general values of fc. Fits to afterglow data us- 
ing fc as a free fitting parameter have yi elded values markedl y 
different from both fc = and fc = 2 ^Starling et al.l 120071 ) . 
although with error bars not excluding either option. The 
scaling coefficients have been obtained by application of our 
post-process code not to a full hydrodynamic simulation but 
to an emulation of this. From the spherical Blandford & Mc- 
Kee (BM) analytical solution for the blast wave for the im- 
pulsive energy injection scenario, snapshots containing the 
state of the fluid at given emission times were constructed 
and stored to provide the input for the post-process code. 

The use of the BM solution provides us with an op- 
portunity to check the results and the consistency of the 
code in an environment where we already have a lot of an- 
alytical control and understanding. The scaling coefficients 
are presented in section [4] They can be used by observers 
to obtain the physical parameters for the blast wave (e.g. 
explosion energy and circumburst density) from the values 
for the peak flux and break frequencies that have been ob- 
tained from fits to the data. Readers interested only in the 
coefficients can skip ahead to this point. The fluxes in the 
transitional regions between the different power law regimes 
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have often been described using heuristic equations that 
smoothly change from one dominant power law to the next. 
The abruptness of this change depends on a sharpness pa- 
rameter s. Using the detailed results from our simulations, 
in section [5] we provide equations for s in terms of two fit 
parameters: the slope of the accelerated particle distribu- 
tion p and the aforementioned k that describes the circum- 
burst density structure. In section [6] we apply our results 
to GRB 970508. We discuss our results in section [7] Some 
cumbersome equations and derivations have been deferred 
to appendices. 



2 DESCRIPTION OF THE 
POST-PROCESSING CODE 

The code takes as input a series of snapshots of relativis- 
tic hydrodynamics configurations on a (in this paper, one- 
dimensional) grid. Although we will treat only the analyt - 
ical Blandford-McKee solution (|Blandford fc McKedll976h 
for the blast wave dynamics put on a grid here, the code 
is written with the intention to interact with the A MRVAC 
adaptive mesh refinement code (jMeliani et alj|2007f ) and will 
read from file the following conserved variables: 
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with 7 the Lorentz factor, p' the proper density, h! the rel- 
ativistic (i.e. including rest mass) enthalpy density, v the 
proper velocity, p' the pressure and c the speed of light. 
From the conserved values we can reconstruct all hydrody- 
namical quantities using the equation of state 



(r Qd - l)e' th , 



(2) 



where F a d the adiabatic index that is kept fixed and e' th the 
thermal energy density. In the entire paper, all comoving 
quantities will be primed. 

The grids represent a spherically symmetric fluid con- 
figuration and all grid cells are assumed to emit a fraction 
of their energy as radiation. This fraction of course has to 
be small enough not to affect the dynamics, since the post- 
processing approach does not allow for feedback. For the 
time being we restrict ourselves to the optically thin case. 

Four ignorance parameters are provided to the code at 
runtime: p, £jv, £e and es, denoting respectively the slope 
of the relativistic particle distribution, the fraction of parti- 
cles accelerated to this relativistic distribution at any given 
time, the fraction of thermal energy that is carried by the 
relativistic electrons and the fraction of thermal energy that 
resides in the (tangled-up) magnetic field. To be precise: the 
fractions €e and es are fractions of e' th , which is strictly 
speaking the sum of the thermal energy of the protons and 
non-accelerated electrons plus the energy of the accelerated 
electrons plus the magnetic field energy. Since we consider 
fully relativistic gases, the adiabatic indices of the electrons 
and protons are both at T a d = 4/3. Also, if the magnetic 
flux enclosed by the surface of any arbitrary fluid element is 
an adiabatic invariant, we find that B 2 oc p 4//3 , which tells 
us that the behaviour of the magnetic energy density B 2 /8iv 
is identical to that of the thermal energy. Or in other words, 
tB retains a constant value away from the shock front. The 
fraction of shock-accelerated particles £jv is often set to one, 
but we have already kept it explicit in our calculations. At 



late times (i.e. when the fluid flow is no longer relativistic) 
£jv has te be lower than unity in order to have enough energy 
per accelerated particle for synchrotron emission. 

In this work we consider synchrotron radiation only. All 
grid cells contain a macroscopic number of radiating parti- 
cles and the radi ation from these p artic le distributions is cal- 
culate d following lSari et al.l l|l998l ) and lRvbicki fc LightmarJ 
(1979), but with two important differences: the transition 
to the lab frame is postponed as long as possible and no 
assumption about the dynamics of the system is used any- 
where as this should be provided by the snapshot files. 

For clarity of presentation we will ignore the effect of 
electron cooling in this section. 

For the emitted power per unit frequency of a typical 
electron we have 
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Here q e denotes the electron charge, m e the electron mass 
(later on we will also encounter the proton mass m p ) and B' 
the local magnetic field strength. The function Q contains 
the shape of the spectrum. It shows the expected limiting 
behaviour: Q(x) oc x 1 ^ 3 for i< 1 and Q(x) oc 2;' 1_p ^ 2 for 
x 1. It incorporates an integration over all pitch angles 
between electron velocities and the local magnetic field and 
an integration over the accelerated particle distribution. We 
use a power law particle distribution with a lower cut-off 
Lorentz factor 7m . Equation @, the critical frequency u' crm 
and the full shape of Q are derived in appendix lAl 

Assuming isotropic radiation in the comoving frame, we 
arrive at 



du' dfi 



4tt du' V ; 



(4) 



per solid angle Q' . 

To get to the received power per unit volume in the 
lab frame, we have to apply the correct beaming factors, 
Doppler shift the frequency and multiply the above result 
for a single particle with the lab frame particle density: 



d 2 P v 
dudQ. 
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with now denoting the cosine of the angle between the 
fluid velocity and the observer (unprimed, so measured in 
the lab frame), ft the fluid velocity in units of c and n the 
number density. 

Finally, the flux the observer receives at a given observer 
time is given by 



du dtt 



(u'(u))(l -Pn)cdAdt e 



(6) 



Here r (, s is the observer distanc^E approximately the same 
for all fluid cells (though the differences in arrival times are 
taken into account). The area A denotes the equidistant sur- 
face. For every emitting time t e a specific intersecting (with 
the radiating volume) surface exists from which radiation 
arrives exactly at t bs- The integration over the emission 
times t e (represented in the different snapshot files) requires 



1 For cosmological distances r j, s denotes the luminosity distance 
and redshift terms (1 + z) need to be inserted in the appropiate 
places in the equations. 
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an extra beaming factor and a factor of c to transform the 
total integral to a volume integral. 

To perform the surface integrals, the post-processing 
code uses a Monte Carlo integration algorithm with both im- 
portance and stratified sampling, using the pseudo-random 
Sobol' sequence^ For the integral over emission times, a 
combination of modified midpoint integration and Richard- 
son extrapolation is used (the latter allowing us to occasion- 
ally skip a snapshot if the desired convergenc e is already 
reache d). All methods are explained in detail in lPress et al.l 
(|l992h . A minor complication is here that not every t e 
probed has a corresponding snapshot file available and inter- 
polation between snapshot files may be needed. The bound- 
aries for surface A are analytically known conic sections and 
depend on the jet opening angle and observer angle. Two 
useful consistency checks are observing a spherical explo- 
sion from different angles and calculating the volume of a 
grid snapshot via integration over different observer times 
while setting the emissivity to one. 

When creating snapshot files directly from the BM so- 
lution we found that sufficient convergence (below the cool- 
ing break) was obtained during the post-processing even for 
modest grid resolutions^ 

For spherical explosions we used jets with an opening 
angle of 180 degrees, which makes no noticeable difference 
for the resulting signal because of relativistic beaming. It 
is worth emphasizing that it is our method that allows for 
the modest grid resolution and keeps calculation time short. 
This is because instead of binning the output from all grid 
cells, it takes an observer time as the starting point and 
then probes the appropriate contributing grid cells only (re- 
solving the structure within the cell by including neighbour- 
ing cells in the interpolation). We have checked our results 
by increasing the accuracy (e.g. larger number of grid cells, 
more snapshot files, smaller step sizes in the integrals etc.) 
and by replacing the Monte Carlo integration routine with 
a nested one- dimensional Bulirsch-Stoer algorithm. These 
consistency checks are in addition to the two mentioned 
earlier. Finally we have checked the grid interpolation and 
snapshot I/O routines by comparing the results of our post- 
process code with those of a code that does not read profiles 
from disc but calculates the BM solution at run time. 



3 THE INCLUSION OF ELECTRON-COOLING 

The code as described so far is purely a post-process code 
that in principle can be applied directly to the output of 
any RHD simulation. If we want to include electron cooling 
however, we can no longer reconstruct the electron energy 
distribution from the conserved quantities alone. In the par- 
ticle distribution function, in addition to the lower bound- 
ary 7 m , we will also have an upper boundary 7m beyond 



2 But if symmetry allows (e.g. the observer is on the jet axis), we 
just do a straightforward Bulirsch-Stoer integration 

3 On the order of 120 base cells with 8 levels of refinement (an 
increase in refinement means a local increase of resolution by a 
factor of two) for a region ^ 10 17 cm to ^ 10 18 cm and a relatively 
small number of snapshots (^ 1000) to go from T ^ 100 down to 
r Ki 2. Unfortunately, the resolution will eventually be dictated 
by that required by RHD simulations, which will be much higher. 
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Figure 1. different possible spectra 



which all electrons have cooled. The time evolutions of both 
the lower cut-off Lorentz factor and the upper cut-off 
Lorentz factor *y' M (that we have tacitly kept at infinity in 
the previous section) of this distribution are no longer dic- 
tated by adiabatic cooling alone but also by radiation losses. 
This implies that when running an RHD simulation we need 
to keep track of at least one extra quantity (at least y' M , al- 
though in practice we will trace both). 

With the introduction of a second critical frequency 
v' cr Mi the equation describing the total emitted power now 
becomes, 



GlP<e> 
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instead of eq. Q. The function Q(xm ,x m ) and 

^cr, m are 

derived and described in appendix [B] For 7m at infinity we 
have Q(0,Xm) — > Q(x m ). 

The particle distribution that lies beneath the deriva- 
tion of this new function Q is no longer a simple power 
law, but drops off sharply for particle Lorentz factors ap- 
proaching the peak value of 7 M . A subtlety worth noting 
here is that the critical frequency v' crM corresponding to 
7m is not the cooling frequency, but a frequency beyond 
which the signal will drop exponentially. Since we put 7m 
at infinity directly behind the shock, we will not directly ob- 
serve v' cr M . The actual cooling frequency is found between 
v'cr.m. and v'cr, My at the point where the shape of the par- 
ticle distribution ceases to be characterized by a power law 
but starts to be characterized by the strong drop towards 
7m- We will discuss the distinction between the cooled and 
uncooled region in appendix iDl 

A consequence of electron cooling is that the amount 
of energy in the shock-accelerated electrons is no longer a 
constant fraction of the thermal energy. e_e now refers to 
the fraction of thermal energy in the shock accelerated elec- 
trons directly behind the shock front instead and the further 
evolution of the available energy is traced via y m and 7m- 



4 SCALING COEFFICIENTS 

Especially for high Lorentz factors, the shape of the spec- 
trum is dominated by the radiation coming from a very thin 
slab right behind the shock front. So we expect the flux to 
scale as 



F <x{p-l)-N t0 f{- 



dfi 



■B'-Q 



(8) 



'(1 — /3/i) 3 7 3 \l^cr,M Vcr,m 

Here N tot is the total number of radiating particles and dfi 
reflects the increasing visible size (due to decrease of beam- 
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ing) of the slab. The two possible spectra that the code can 
generate are shown in fig.[T] where we used the labelling from 
iGranot fc~S ari (2002) to distinguish the different power law 
regimes. In tables {TJ and ((2)1 we give the expressions for 
the absolute scalings in the different regimes D, E, F, G, H 
and the critical frequencies. Scaling co efficients aside, these 
equat ions are similar to those given in IVan der Horst et al.l 
( 2008). The flux in regime D is denoted by Fd, the critical 
peak frequency in spectrum 1 is denoted by v m> \, the critical 
cooling frequency in spectrum 1 by v Ci \ and so on. 

The equations in the tables introduce a number of sym- 
bols that need an explanation. The cosmic redshift is given 
by z, while the luminosity distance r b s , 2S is measured in 
units of 10 28 cm. E52 is the explosion energy E in units 
of 10 52 erg. The observer time in days is denoted by t bs.d- 
The characteristic distance Ro we put at 10 17 cm and po and 
no are related via the proton mass: po = m p no. The scal- 
ing coefficients Cd, Ce etc. contain a number of numerical 
constants (determined by fitting to output from our code) 
and some explicit dependencies on k and p and are further 
explained in appendix [Cl 

Before the cooling break the scaling behaviour is dic- 
tated by the asymptotic behaviour of Qiv' lv' crm ). The 
steepening of the spectrum beyond the cooling breaks and 
the corresponding changes in the scaling behaviour are due 
to the fact that beyond the cooling break frequency the re- 
gion behind the shock that still significantly contributes to 
the total flux (i.e. the hot region) becomes noticably smaller 
than the shock width. The changes in the scalings reflects 
the change in the size of region. The hot region is discussed 
separately in appendix [Dl 



5 SHARPNESS OF BROKEN POWER LAW 

In simple power law model fits, the gradual transition be- 
tween regimes is often handled by a free parameter, the 
sharpness factor s. In more detailed calculations like those 
done here the gradual transitions are included automatically 
and we can use this to provide the correct dependence of s 
on p and k. This eliminates s as a free parameter, simplifying 
the fit to the data and allowing the shape of the transition 
to help determine whether a particular model fits the data 
or not. 

For spectrum 1, we use the following equation to de- 
scribe the flux density near the peak break v m ,i- 



few percent: 



F{v) = F m ,i- 



i.i(i-p) 



v 



+ 



V 



(9) 

where F m ,i denotes the flux at the critical frequency v m ,i for 
infinite sharpness s m ,i (i.e. the meeting point of the asymp- 
totic power laws) . When we switch off cooling in our simula- 
tion, we can determine s m ,i from fitting against the result- 
ing spectrum while keeping the other parameters in equation 
((9| fixed. The sharpness is a function mainly of p and to a 
lesser extent of k and the other simulation input parameters. 
Rather than attempting to include all secondary dependen- 
cies when formulating a description for s m ,i, we find that 
the following approximation for s m ,i is always valid up to a 



s m ,i = 2.2 - 0.52p. 



(10) 



When we switch on electron cooling, the flux is best approx- 
imated by 
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If we fit this function against simulation output using s e ,i 
as a fitting parameter we find that the results are described 
(up to a few percent) by 



s c ,i = 1.6 - 0.38p - 0.16/c + 0.078pfc. 



(12) 



A simultaneous fit using both s m ,i and s c ,i yields the same 
results. 

For spectrum 5 the order of the breaks is reversed and 
the smooth power law for both breaks is given by 

1 



F(u) 
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(13) 



where F C) g denotes the peak flux for infinite sharpness s Ci5 
and the prescriptions for the sharpness are 



and 



s c , 5 = 0.66 - 0.16fc, 



s m , 5 = 3.7 - 0.94p + 3.64fc - 1.16pk. 



(14) 



(15) 



Once again valid up to a few percent. Given their ac- 
curacies, all sharpne ss prescriptions are consistent with 
IGranot fe Saril (|2002l ). 



6 APPLICATION TO GRB 970508 

Various authors have used flux scaling equations to derive 
the physical propert i es of GRB 970508 from afterglow data 
ilGalama et al J 1 19991: IGranot fc Saril 120021 : lYost et ail 120031 : 
IVan der Horst et al. 2008). This provides us with a context 
to illustrate the scaling laws derived in section|4] We will use 
the fit parameters obtaine d from broadband modeling by 
IVan der Horst et ail (|2008h . They have fit simultaneously in 
time and frequency while keeping k as a fitting parameter. 
Because the only model dependencies that have been intro- 
duced by this approach are the scalings of t and v (and no 
scaling coefficients), their fit results are still fully consistent 
with our flux equations. Using the cosmology Qm = 0.27, 
S7a = 0.73 and Hubble parameter Ho = 71 k m s -1 Mpc" 1 , 
they have r b s ,28 = 1.635 and z = 0.835 jMetzger et ail 
1 1997ft . leading, at t obsA = 23.3 days, to u cA = 9.21 • 10 13 
Hz, Um,i = 4.26 ■ 10 10 Hz, F m ,i = 0.756 mjy, p = 2.22 and 
k = 0.0307. 



Gamma-Ray Burst afterglow scaling coefficients for general density profile 5 

Table 1. Flux Scalings for the different regimes (see tables ICT1 and IC2l for C D ,C E etc.) 
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Table 2. Critical frequencies for the different regimes (see tables ICll and IC2l for Cb,Ce etc.) 
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Both IVan der Horst et all (|2008h and Gal ama et al.l 
( 1999) take for the hydrogen mass fraction of the circum- 
burst medium X = 0.7, which in our flux equations is 
mathematically equivalent (though conceptually different) 
to setting £jv = (1 + X)/2 = 0.85. Unfortunately this still 
leaves us with four variables to determine (en, £b, -E52, no) 
and only three constraints (peak flux, cooling and peak fre- 
quency). From a t heoretical study o f the microstructure of 
collisionless shocks iMedve dcv (2006) arrives at the following 
constraint: 



(16) 



We include this constraint to have a closed set of equations. 

For the values quoted above we obtain: £52 = 0.155, 
n = 1.28, e s = 0.1057, e E = 0.325. In figures [2] and |3] we 
have plotted a comparison between the spectrum generated 
by using these values as input parameters for the BM solu- 
tion and the spectrum as it is represented by apply i ng the 
results of the broadband fit of IVan der Horst et all (|2008h 
for the values of the critical frequencies and the peak flux 
to equation (111[) . 

Our scaling coefficients were fixed for arbitrary k and 
for comparison we also give results for k — and k = 2. 
The ISM case is virtually identical to k — 0.0307 and yields: 
£ 52 = 0.155, n = 1.23, e B = 0.106, e E = 0.325. The stellar 
wind case yields: £52 = 0.161, no = 6.45, en = 0.0957, 



GRB970508 spectrum after 23.3 days 
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Figure 2. A comparison of a representation of the data, using 
the values of Van Der Horst for the critical frequencies and peak 
flux, combined with our equations for the transition sharpnesses, 
at 23.3 days, with a reproduced spectrum from a BM blast wave 
simulation. 



e E — 0.309. The quantity no (the particle number density 
at the characteristic distance 1 ■ 10 17 cm) is affected most. 
Also for comparison we give some of the values ob- 
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simulated light curves GRB970508 
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Figure 3. Simulated light curves generated from post-processing 
of simulations using the BM solution with the input parameters 
we have derived: E 52 = 0.155, n = 1.28, e B = 0.105, e E = 0.325. 
For the X-rays we used v = 1 ■ 10 18 Hz and for the R band 
v = 4.3 ■ 10 14 Hz. 



tained by other authors. iGalama et al.l (|l999h obtain for the 
ISM case: E 52 = 3.5 , n = 0.03, e fl = 0.09, e E = 0.12. 
iGranot fc Saril (|2002T l obtain for the ISM case: E 52 = 0.12, 
n = 22, e B = 0-012, e E = 0.5 7. Both use p = 2.2. Fi- 
nally IVan der Horst et all (120081 ) obtain for k = 0.0307: 
E &2 = 0.435, n = 0.0057, e B = 0.103, e E = 0.105. 

The large differences between the various results illus- 
trate the importance of using the correct scaling coefficients 
to derive physical parameters of GRBs and provide a strong 
motivation for this work. Because the error ones in particu- 
lar is rather large for the quoted authors, who have used the 
self absorption critical frequ ency to pr o vide a fourth con- 
straint, the constraint from iMedvedevi (|2006l ) can not be 
rejected based on their fit results. The extension of our code 
to include self-absorption will yield an alternative and can 
be used to further study the applicability of Medvedev's 
constraint. 



7 SUMMARY AND DISCUSSION 

In this paper we have introduced an approach to reconstruct 
light curves and spectra from hydrodynamic simulations. 
The central idea is that we do not start from simulation 
snapshots and bin the output of each grid cell, but that for 
representative snapshots we integrate over the intersecting 
surface that contains all points where radiation is gener- 
ated that is due to arrive at a given observer time and fre- 
quency. When performing these integrations we interpolate 
within and between grid cells. While in the context of this 
paper we have used only snapshots that contain mimicked 
RHD output using the BM solution, first results using real 
simulations have been obtained and will be discussed in a 
later paper. An important thing to note here is that, even 
though the post-process code only required a very modest 
resol ution, the underlyin g hydro dynamics code usually does 
not. iMeliani fc Keppensj (|2007t) used 1200 base level cells 
and 15 refinement levels to simulate the evolution of the 



blast wave (earlier, when they were putting their code to 
the test they even used 30,000 base level cells at one point, 
see IMeliani et aT1l2007l ) . This means that, in general, paral- 
lel computer systems are required to run these simulations, 
something for which the RHD code that our post-process 
code interacts with (AMRVAC) was explicitly designed. 

In our code we included synchrotron radiation and elec- 
tron cooling. We use a parametrisation of the accelerated 
particle distribution in terms of and 7,„. Thermal ra- 
diation from the particles not accelerated to a power law 
distribution can be included in a straightforward manner. 
The code can also be extended to include self-absorption 
and since the outgoing synchrotron radiation from a grid 
cell is independent of the incoming radiation, this can be 
done without expanding to a full radiative transfer code in- 
cluding scattering. Effectively, all that is needed is to post- 
pone the integration over the intersecting surfaces until after 
the integration over emission times, while in the meantime 
diminishing the output from earlier surfaces according to 
the column densities in the lines of sight, which amounts to 
solving linear transport equations only. 

As a consistency check and a first application of the 
code we calculated the scaling coefficients of the flux scaling 
equations for GRB afterglow spectra for arbitrary values of 
k with unprecedented accuracy. These results can be used to 
obtain the physical parameters of the burst from fits to af- 
terglow data. For the ISM and stellar wind sc enario's the re- 
sults h ave been checked against the results of lGranot fc Saril 
(2002) and are found to be fully consistent. The motivation 
for the choice of arbitrary k is tha t various authors have now 
used k as a fittin g parameter (e.g. lVan der Horst et al ] |2008l : 
I Yost et al]|2003l ). Values of k other than or 2 reflect the 
structure of a circumburst medium altered by shock interac- 
tions or more complicated stellar wind structures. We have 
used GRB970508 to illustrate the effect of using our scaling 
coefficients to deduce the physical properti es of a GRB. Here 
we have used an additional constraint by IMedvedevi (|2006l ) 
to obtain a closed set of equations in the absence of a full 
description for the self-absorption. 
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APPENDIX A: DERIVATION OF EMITTED POWER PER ELECTRON 

For each electron Lorentz factor -y e we define two critical frequencies v' cr t , a and v C r,e- 

3 ,2 1e.B' i /,,, 

= - — 7 e sma = v cre sma., (Al) 



where q e denotes the electron charge, m e the electron mass and a the pitch angle between field and velocity. It is around (but 
not exactly at) these values that the spectrum peaks and we will find them useful as inte gration variables later on. 
The power per unit frequency emitted by an electron is (IRvbicki fc LightmanI (I1979T )): 

where 

F(x) = x Ks(Od£, (A3) 

J X 

with K5 a modified Bessel function of fractional order. F(x) behaves as follows in the limits of small and large x: 

^^©""O-^'^H' i<<i - <A4) 

w s /f i/ 2 _„ A 55 1 10151 1 \ , 

F(x)^\ -x ' e H , *>1, (A5) 



2 V 72 1 10368 

where T(x) is the gamma function of argument x. 

For the mean power averaged over all pitch angles while assuming an isotropic pitch angle distribution we obtain: 

where 

V(x) = - (smafFi-^- ) da. (A7) 
2 J sm a 

In the limit of small and large x, P(x) behaves as follows: 

V{x) ~ 22/3 • ^ • ^ x 1 ' 3 -I=x- 5 • ^ ^ x 7 / 3 , x«l, (A8) 
9r(^) V3 48-2Vsr(^) v ' 

- | e -", x » 1. (A9) 

The effective lower cut-off Lorentz factor of a collection of electrons 7^ can be expressed in terms of local fluid quantities. 
The integrated power law particle distribution C , y e ~ p d^ (C is a constant of proportionality) must yield the total number 
density of particles: 

jT C( 1 ' e y P d< = Civn' - C = -^T^n'. ( A1 °) 
Similarly the integrated particle energies must yield the total energy: 



(All) 



Combining these equations and dropping the rest mass term in the energy equation (it will be negligible for relativistic 
electrons), we obtain 



2 -p\ fe E e' th 1 



1—pJ \£n n' rn e c 2 



(A12) 



If we integrate ()A6|) over the particle distribution and divide the result by the total electron density, we obtain the emitted 
power per ensemble electrorQ 

^(V) = *=± ■ ■ Q (-^) . (A13) 

dv> y > 2 m e c 2 ^\v> ) K ' 



4 an ensemble electron contribution is therefore constructed as the total of all electron contributions divided by the number of electrons. 
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Figure Al. F(x), "P{x) and Q(x) (for p = 2.2 and p = 2.8 ): single electron, single angle-averaged electron and ensemble electron 
respectively. 



Here v' cr m denotes the resulting value of v' cr t , when we substitute "/' m for y' e in equation (|A1|) . It surfaces when we switch 
integration variables from *y' e to v' cr e . The auxiliary function Q is defined as 



Q(x) 



i-p 

: x 2 



p-3 

v 2 ~P(y)<iy- 



In the limit of small and large x, Q(x) behaves as follows: 

2 5 ^r(|) 1/3 

Q(X) ^ — TT 2 — X 



5(3p - 1) 



-x + 



3\/3^r(| 



V3( P +1) 2V3(88 + 24p) 



6/ x 7/3 , 



x < 1, 



(A14) 



(A15) 



x > 1. 



(A16) 



, _r(| + f) 2— ^.p 19, w p 1. i^p Tre x 

o(x)^Jn—4 — 14 • r(- h — )r(- )-k 2 , 

^ ; v r(| + f) p+1 v 4 12 ; ^4 12 ; 2 a; 

In practice, the computer code uses lookup tables for F(x), V{x) and Q(x). The three functions have been plotted in figure 
[[AT) (Q for both p = 2.2 and p = 2.8), allowing for comparison between the spectra from a single electron, an angle-averaged 
electron and an ensemble electron. 



APPENDIX B: EMITTED POWER WITH ELECTRON COOLING 



If the only processes that are of importance are synchrotron emission and adiabatic cooling, the evolution of the Lorentz 
factor of a single electron is described by 



d 7e = (TT(B') 2 2 | 7e dri 

dt' Q-Km e c ^ e 3n' dt' ' 



(Bl) 



where or denotes the Thomson cross section. In lGranot fc Sar i (2002) this differential equation is applied to the BM solution 
by expressing it in terms of the self-similar variable and solving it analytically. In our case we can use eq. ()A12|) to establish 
7^ directly behind the shock front and initially put rf M , the upper cut-off Lorentz factor due to cooling, at a sufficiently large 
value (instead of infinity). Sufficiently large for example can be taken such that 



/•oo 
v 7 



«S e, 



(B2) 



with e some tolerance for the error in the energy. The real will quickly catch up with the approximated 7^-, as can be 
seen from equation (|B1|) . 

The analytical solution for the particle distribution in the BM case is given by 



N e (<) = Cj'~ p ■ (1 



7e y-2 



7m 



(B3) 



10 H.J. van Eerten and R.A.M.J. Wijers 



where the factor C now stands for 

C = (p - l)en' 7 r 1 • (1 - ^) 1 - p . (B4) 

7m 

We take this to hold for the output of real RHD simulations as well, so that we have an approximate parametrisation for the 
part icle distribution in any g rid cell in terms of y' m and j'm alone. A more complete treatment of the particle distribution 
fe.g. lPe'er fc WaxmarJ ((2005) ) would effectively introduce an additional dimension to the RHD simulation and slow down the 
calculations accordingly. 

Via reasoning completely analogous to the non-cooling case (where we use eq. (|A1|) with j' M instead of j' e to obtain v' cr M ) 
we arrive at an auxiliary function Q given by 

Q(vM, ym ) = y y ■ (i - (^i) 1 / 2 ) 1 - . r y ^.(i- { y^)^r 2 ny) * y , (b 5 ) 

y™ j VM y 



ocurrmg m 



diy M _szi.^w. a| - - , (B6) 



cr,Al "cr,Tn 



Since this is a function of two variables instead of one, its limiting behaviour is more complicated. If yM <C 1, Q(jjM ,Vm) 
approximately reduces to 

Q{y M ,y m ) X {l-(^L)^f-^ -Q(y m ), J/M«l, (B7) 

y™ 



which can be obtained by approximating yu by zero in the integration limits and integrand of equation (|B5[I . If ym/yu —* 1, 
the approximate result is 

Q(yM,ym) oc V(y m ), ^ - 1, (B8) 
Vm 



which follows from approximating the integral from ()B5[) by its value at y m times the integration domain. If y m < 1 as well, 
we can use the first term of the lower limit series expansion for V in the integral and solve it to refine the approximate result 
to 

Q(y M ,y m )(xP(y m )/(p-l), j/ m «l. (B9) 

VM 

On the other hand, if y m /yM S> 1, we can approximate the result in terms of Q(yM,y' m ) for a smaller value y' m (i.e. the 
last tabulated value): 

Q(yM,y m )^Q(yM,y' m )-(^P-) * + Q(ym) - Q{y' m ) i^f] , ^ » 1, (BIO) 
\y m J \y m J vm 

which further reduces to 

Q(yM,y m )~Q(y M ,y m ).{-) --_ + -_. , - » l, (Bll) 

for sufficiently high values of y' m and y m . 

Finally, for yM > 1 we find from fitting to tabulated values that Q(yM ,y m ) is best described by 

i-p 



Q(y M ,y m )oc[^) • (1 - (i^)^)^ . & -vm . y v-^ yM » L (B12) 
\yM J ym 

In practice the code uses a two-dimensional table with numerically calculated values in addition to the analytical expressions 
above. The contribution from the region where yM S> 1 is effectively zero due to the exponential term e~ VM . 



APPENDIX C: DERIVATION OF SCALING COEFFICIENTS 

We summarize the equations for the scaling coefficients in table (1C1[) . Aside from some explicit dependencies on p and k these 
equations also contain truly numerical constants with values that have been determined by fitting to output of our code. For 
example Co{p,k) contains the constants Cdo, Cok and Cokk (with C% h denoting Cok to the power k etc.). Their numerical 
values are listed in table ()C2|I . Instead of incorporating these numerical constants in the total flux formula as we have done 
here we could also have used a fitting polynomial, but this approach more closely reflects the k and p dependencies. 

The first term (p — 1) in these equations is also the first term in eq. (|8]). From the contribution of Ntot we obtain a 
contribution 1/(3 — k) via 

^-^CM-kY (C1) 
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Table CI. scaling coefficients 



C D = (p - 1) ■ 
•(17 - 4fc 


< *~i /~ih /~ik 2 

10— 4k 
)3(4-fc) . (4-fc) 


,V(4-fc) 1 /p-2\~ 2/3 
' 3-fc Vp-1/ 

2-fc 

4=1 .Q_ 


Ce = (p - 1) • 
■(17 - 4fc 


^ ^ ^ 2 y/( 4 - fe > i 

y^EO^Ek^Ekk ) ' 3 — fc 

-efc+14 2-3fc 
) 3(4- fe ) . (4 _ fc) 3(4- fc , .Q cQoi 


C* F = (p - 1) • 
■(17 - 4fc 


c FO c Fk c Fkk ) ■ 3 _ k 

) 3 / 4 ■ (4 - fc)- 1 / 4 ■ Q cooi 


C G s (p - 1) • 
■(17 - 4fc 


.-.fc ,-.k 2 ^-pfc /^pk 2 n p 2 r <p 2 k n p 2 k 2 \ 
^GO<~'Gk L 'Gkk L 'G P L 'Gpk L/ Gpkk u Gpp L/ Gppk L/ Gppkk J 

-kp-5k + 4p+12 3fcp-5fc-12p+12 

) 4(4- fe ) -(4-fc) 4 < 4 - fc > -Q+ 


l/(4-fc) 1 , 
3-fe \ 


'p-2\ 




Cff EE (p - 1) . 

(17 - Ak) 


f (~i (~ik fyk 2 (~<P rypk (~,pk 2 f-tp 2 (~*p 2 k {~*p 2 k 2 
v <- \ffO t -' Hk>- Hkk'~ Hp'~ \Hpfc L HkkS- "Hpp~ HppW~ Hppkk 

p+2 2-3p 

— ■ (4-fc) - a - -Q+ 


l/(4-fc) 1 

' 3 - fc 




r 



Table C2. Constants setting scale of flux 



D G H F E 



5.12 10" 17 2.78 10" 31 5.68 • 10" 1 1.16 • 10 30 2.95 ■ 1CT 16 

fc 1.18 ■ 10 4 4.54 10 7 6.94 • 10" 1 1.36 • 10~ 8 2.04 ■ 10 4 

fcfc 9.01 ■ lCT 1 8.95 10" 1 9.27 • 10" 1 1.01 9.41 ■ 10" 1 

p 2.25 ■ 10 32 5.40 • 10 30 

pk 7.27 ■ 10~ 9 1.65 • 10" 8 

pfcfc 9.41 ■ 10" 1 1.06 

pp 1.77 2.99 

ppk 8.07 ■ 10" 1 7.01 ■ 10" 1 

ppfcfc 1.03 1.01 



The o rigin of the combination (p — 2)/(p — 1) can be traced to equation (|A12|) in appendix [Al of this paper (jGranot fc Sari 
(2002)) actually absorb it into ce). The term (17 — Ak) is linked to the energy £52 and the two will always occur with the same 
power as can be seen from comparing tables IC II and [1] It enters our calculations via equation (69) from iBlandford fc McKeel 
l|l976l) . The term (4 — fc) is likewise linked to the observer time t t s ,d- The extra term is a result from the transition from 
emission time in the grid lab frame to observer time. For the shock front the two are related via 

t e = (2(4-fe)U s ) [ 8npoc5 -k R k) • (C2) 

The final terms are different for the different power law regimes. They are contributed by the leading order terms of the 
various approximations of Q. Q+ is given by (see eqn. (|A16|l ): 



r(5 + E)r(E + ±2)r( E - ±) 

q + — L V4 ~ 4>^ V4 ~ 12 ^ U 12' 

For low uncooled frequencies we have 



r(I + f)(p + i) 



(C3) 



(C4) 



3p - 1 

as can be seen from equation (|A15|I . When cooling plays a role we find that equation (|B9[) provides us with 

Qcooi=— l —- (C5) 

p- 1 

Note that the effect of Q coo i in Ce and Cf is to cancel out the first (p — 1) term -we only kept both terms for clarity of 
presentation. 
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APPENDIX D: THE HOT REGION 

For any given observer time and observer frequency there is a region behind the shock front where the emitting electrons 
have not yet cooled below the observer frequency. Although, when we set 7m initially at infinity, the size of this region never 
becomes zero, it can become very small, even when compared to the analytical error on the BM solution. The size of the hot 
region also determines the slope of the spectrum beyond the cooling break. We calculate its properties below. 

The BM solution is obtained by a change of variables from t and r to \ and 1 /F 2 , where the fact that the latter becomes 
very small is continually used to simplify the dynamic equations using first order approximations. The % coordinate of a fluid 
element is given by 

X =[l + 2(4-fc)r 2 ](l-^), (Dl) 

which is 1 at the shock front and increases roughly one order in magnitude until the back of the shock. 

The radiation received at a given observer time is obtained by integrating over equidistant surfaces that have a one-on-one 
correspondence to emission times. To obtain an order of magnitude estimate for the size of the hot region we look solely at 
the jet axis, where each emission time and hence each equidistant surface corresponds to a given position Xi via 

X (r,t) = X (c(t e - t obs ),t e ) « tfl . 2(4 - fc)B 2 . (D2) 

We define the boundary of the hot region Xhot at the point where u' cr M = v' (i.e. when the second argument of Q(—r — , —r- — ) 

is equal to one). The critical frequency z^ r ,Af is related to via the u s ual re lation (see eqn. lAlj) . and an expression for -y' M 
in terms of the self-similar parameter % can be found in iGranot fc Saril l|2002l ): 

i ( \ _ 2(19 - 2k)nm e c~/ 1 , . 

7mW- a T BHl x (19-2fc)/3(4-fc) _ I ■ { - LKi ) 



Using the above we can find an expression for Xhot -or equivalently t e: hoti since the two are related via eqn. (|D2[) . To first 
order in Xhot — 1 we find 

f 27g e m e (4-fc) 2 \ 1/2 / E(17 - 4k)t obs 2(4 - k) \T0h 

Xh0t ~ \vall2^cHT P l /2 Rf /2 ) V 8tv P0 R^ J [ > 

1 / 27g e m e (4-fc) 2 \ f E(17 - 4k)t obs 2(4 - k) \^=V 



where t e j r0 nt is the emission time of the shock front (see equation 102ft . 

From this we can draw a number of conclusions. The size of the hot region is dependent on the observer frequency via 
j,- 1 / 2 p or observer frequencies beyond the cooling break, it is effectively this region alone that contributes to the observed 
flux, since the contribution from the cool region drops exponentially (see equation IB 12fl . A steepening of the spectral slope 
by -1/2 is therefore expected: (1 — p)/2 — + —p/2. This results from multiplying the pre-cooling break flux by the fraction of 
the total emitting region that consists of the hot region -which is given by 



Ujront - t e ,hot J / 27g e m e (4- fc) 2 \ 1/2 /£(17-4fc)t o6s 2(4-fc)\4ni 



ie,/ r0 ni A—k \^l28V^^e 3 J 2 p 3 /2 R 3 k/2 J V 8np R%c 5 - k 



(D6) 



Note that from this equation all post-cooling break scalings (e.g. e_g, £52 etc.) can be derived by multiplying with the relevant 
pre-cooling flux. 

Another important issue is that the size of the hot region can become smaller than the analytical error inherent in the BM 
solution, which cuts off beyond 1/F 2 . This happens at late times, when T has dropped significantly. This puts a practical limit 
on a direct numerical implementation of the BM solution in our radiation code, ironically not due to numerical limitations 
of the code but because of the upper limit on the accuracy of the analytical solution that we have used to generate our grid 
filefl On can however still extrapolate the heuristic description of the spectra and light curves that we have obtained for 
arbitrary fc to late times -this is completely consistent with the canonical approach to light curve and spectrum modelling. 

This paper has been typeset from a TgX/ F/IpX file prepared by the author. 



5 In general, when post-processing grid files from simulations the issue does not occur because we use the AMR structure of the grid 
to set the local integration accuracy. If the hot region becomes very small, then this will be dealt with at the earlier stage of the RHD 
simulation. Also, when directly integrating the flux equations for the BM solution by first expressing everything in terms of the self-similar 
coordinate and sticking to that frame, the issue is largely avoided as well. 



